/**
 * @file LMM.h
 * @brief The class for LMM.
 * @author XDDDD
 * @version 
 * @date 2021-04-28
 */

#ifndef __PAISLEYPARK__LMM_H__
#define __PAISLEYPARK__LMM_H__

#include "TimeIntergrator.h"
#include "LMM_Info.h"
#include <cstring>
#include <fstream>
#include "RK_Info.h"

template <std::size_t DIM>
class LMM : public TimeIntergrator<DIM> {
public:
	/**
	 * @brief Set the coefficients of the lmm.
	 *
	 * @param _file_in Load data from here.
	 * @param _p The given order of p.
	 */
	void set_cof(std::string _file_in, const std::size_t _p);
	/**
	 * @brief Get the steps of the lmm.
	 *
	 * @return number of step of the lmm.
	 */
	const std::size_t get_s() const;
	/**
	 * @brief Get the order of the lmm.
	 *
	 * @return order of the lmm.
	 */
	const std::size_t get_p() const;
	/**
	 * @brief Get the number of initial points for each step.
	 *
	 * @return The number of initial points.
	 */
	const std::size_t get_num_init() const;
	/**
	 * @brief Get alpha_i
	 *
	 * @param _i The given index.
	 *
	 * @return alpha_i
	 */
	const double get_alpha_i(const std::size_t _i) const;
	/**
	 * @brief Get beta_i
	 *
	 * @param _i The given index.
	 *
	 * @return beta_i
	 */
	const double get_beta_i(const std::size_t _i) const;
	/**
	 * @brief Reset the information of lmm.
	 *
	 * @param _a The given alpha list.
	 * @param _b The given beta list.
	 * @param _s The given step number.
	 * @param _p The given order.
	 */
	void set_info(const double* _a, const double* _b, const std::size_t _s, const std::size_t _p);
	/**
	 * @brief Calculate the initial points from one given point by Runge-Kutta of order p.
	 *
	 * @param _Pnt The given vector of points.
	 * @param _stepl The length of step.
	 * @param _num The number of initial points needed.
	 * @param _f The given vector functions.
	 */
	void starting_method(std::vector<Point<DIM> >& _Pnt, const double _stepl,
						 const std::size_t _num, VectorFunction<DIM, DIM>& _f) const;
protected:
	/**
	 * @brief Store the coefficients of lmm of order p.
	 */
	LMM_Info __LMM_Info;
	/**
	 * @brief Store the coefficients of Runge-Kutta of order p.
	 */
	RK_Info __RK_Info;
	/**
	 * @brief The number of initial points needed for this lmm.
	 */
	std::size_t __num_init;
};

template <std::size_t DIM>
void LMM<DIM>::set_cof(std::string _file_in, const std::size_t _p) {
	std::size_t order = 0;
	std::size_t s = 0;
	std::size_t num_init = 0;
	bool is_found = false;
	std::ifstream data_in(_file_in);
	std::string tmp;
	std::getline(data_in, tmp);
	while(!is_found && !data_in.eof()) {
		data_in >> s >> order;
		if(order < _p) {
			std::getline(data_in, tmp, '%');
		}
		else {
			data_in >> num_init;
			if(num_init < s) {
				std::cerr << "Do not have enough number of initial points." << std::endl;
				std::exit(-1);
			}
			__num_init = num_init;
			data_in >> tmp;
			double a[s];
			double b[s];
			double x1,x2;
			for(int i = 0; i != s + 1; i++) {
				data_in >> x1 >> x2;
				a[i] = x1/x2;
			}
			data_in >> tmp;
			for(int i = 0; i != s + 1; i++) {
				data_in >> x1 >> x2;
				b[i] = x1/x2;
			}
			is_found = true;
			this -> set_info(a, b, s, order);
		}
	}
	if(!is_found) {
		std::cerr << "Cannot found such lmm that satisfies the required accuracy." << std::endl;
		std::exit(-1);
	}
	data_in.close();

	std::ifstream rk_in("RK_ButcherTableau");
	is_found = false;
	int pole = 0;
	int step_num;
	while(!is_found && !rk_in.eof()) {
		double x1, x2;
		rk_in >> tmp >> order;
		if(order < _p) {
			std::getline(rk_in, tmp, '%');
		}
		else {
			is_found = true;
			rk_in >> tmp >> step_num;
			double c[step_num];
			double a[(step_num - 1)*step_num/2];
			double b[step_num];
			int pole = 0;
			for(int i = 0; i != step_num; i++) {
				rk_in >> x1 >> x2;
				c[i] = x1/x2;
				for(int j = 0; j!= i; j++) {
					rk_in >> x1 >> x2;
					a[pole] = x1/x2; 
					pole++;
				}
			}
			for(int i = 0; i != step_num; i++) {
				rk_in >> x1 >> x2;
				b[i] = x1/x2;
			}
			__RK_Info.set_info(a, b, c, step_num);
		}
	}
	rk_in.close();
	if(!is_found) {
		std::cerr << "Can not find the RK meeting the given order." << std::endl;
		std::exit(-1);
	}
};


template <std::size_t DIM>
const std::size_t LMM<DIM>::get_s() const {
	int s = __LMM_Info.get_s();
	return s;
};

template <std::size_t DIM>
const std::size_t LMM<DIM>::get_p() const {
	int p = __LMM_Info.get_p();
	return p;
};

template <std::size_t DIM>
const std::size_t LMM<DIM>::get_num_init() const {
	return __num_init;
};

template <std::size_t DIM>
const double LMM<DIM>::get_alpha_i(const std::size_t _i) const {
	double x = __LMM_Info.get_alpha_i(_i);
	return x;
};

template <std::size_t DIM>
const double LMM<DIM>::get_beta_i(const std::size_t _i) const {
	double x = __LMM_Info.get_beta_i(_i);
	return x;
};

template <std::size_t DIM>
void LMM<DIM>::set_info(const double* _a, const double* _b, const std::size_t _s, const std::size_t _p) {
	__LMM_Info.set(_a, _b, _s, _p);
};

template <std::size_t DIM>
void LMM<DIM>::starting_method(std::vector<Point<DIM> >& _Pnt, const double _stepl,
							   const std::size_t _num, VectorFunction<DIM, DIM>& _f) const {
	if(_Pnt.size() != 1) {
		std::cerr << "The number of starting point is not 1!" << std::endl;
		std::exit(-1);
	}
	Point<DIM> Un = _Pnt[0];
	int pole;	
	Point<DIM> y[__RK_Info.get_s()];

	for(int i = 1; i != _num; i++) {
		pole = 0;
		for(int j = 0; j != __RK_Info.get_s(); j++) {
			y[j] = Un;
			for(int k = 0; k != j; k++) {
				y[j] = y[j] + __RK_Info.get_a_i(pole)*_stepl*y[k];
				pole++;
			}
			y[j].set_time(Un.get_time() + _stepl*__RK_Info.get_c_i(j));
			y[j] = _f(y[j]);
		}
		for(int j = 0; j != __RK_Info.get_s(); j++) {
			Un = Un + _stepl*__RK_Info.get_b_i(j)*y[j];
		}
		Un.set_time(Un.get_time() + _stepl);
		_Pnt.push_back(Un);
	}
}
#else
//Do nothing.
#endif
